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Abstract 
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I. INTRODUCTION 



Numerical relativity has undergone a rapid development in the past few years. After the break- 
through of |Q1 HHH], stable longterm simulations of binary black hole (BBH) systems are common 
practice, besides waveform modelling, to study the close-to-mer ger spin p recession or to 

model the final spin [fi 0, S, @, GlJ of BBH inspirals (TIL QJ, QJ, Q3, [3 03]. Recently exten- 
sive investigations have been done concerning the formation process and spin evolution of black 
holes with accretion disks Inl [lgjll appe aring in fully relativistic simulations of binary neutron 



24], rotating neutron star collapse 11251 1261 1271. 12811 and 
3. 



stars [|19l.l20l.l2 111, mixed binaries |2J 
rotating supermassiv star collapse |2i 

In these cases accurate numerical techniques to extract the spin of a BH in a gauge invariant 
manner are required. It is common to obtain a rough approximation of the spin through the quasi- 
normal mode oscillation extracted from the gravitational waveform after merger within black hole 
perturbation theory. Another approximation scheme is to integrate the radiated angular momentum 
contained in the gravitational radiation at 'large' coordinate spheres to draw conclusions about the 
remaining spin of the system given the initial data. 

Other methods, as discussed in this paper, use the gauge invariant notation of an apparent 
horizon (AH) or in more general terms a marginally outer trapped surface (MOTS) which can be 
located on the spatial slices of the simulation. There gauge invariant spin and mass can be defined, 
if an axial Killing vector field (KVF) $ a is present, as in the case of Kerr. But opposed to the 
stationary case, the spacetime outside the horizon can be dynamical without spoiling the gauge 
invariance of these quantities !132ll33L l34l. 13511 . The invariant quasi-local spin J[$ J ] is given by the 
surface integral (Brown- York form) 



J[$J] : = -— (f) &K ijS l dA. 



57T 



(1.1) 



where dA is the 2D area element, the extrinsic curvature of the Cauchy slice and s % is the 
outward-pointing surface normal on the MOTS denoted by S. In order to obtain the 2D Killing 
equation has to be solved; if the axisymmetry is perturbed approximate KVFs (aKVFs) have to 
be computed [36[ 37, 38], for applications in BBH simulations see [11, 16]. Sometimes, due to 
computational reasons, the effort of finding a KVF or aKVF is not done and coordinate vector 
fields are instead used to estimate J[$ J ] « ^[^ v ]> see e.g. |0, H^. Another common set of 
methods to determine the spin uses properties of the Kerr solution at the horizon, such as the 
proper length of the 'equatorial' circumference rf40Tl or the extrema of the scalar 2-curvature [16]. 

In this paper we present a new, comparatively easy to implement algorithm, which is based on a 
multipole decomposition of the rotational Weyl scalar Im\l/2 114 ill in the framework of the isolated 
and dynamical horizon formalism |33|,[3l],[35f|; for reviews see e.g. IsHEIH]. The dipole term 
reads 

Ji = -\H^4- 2 Y w ( X )dA, (1.2) 



12vr 4vr 



where A is the horizon area, (x, 4>) an invariant coordinate system |41|] 'tied' to the axisymmetry, 
such that Ji and J[$ J ] are identical, and Y 10 (x) is the spherical harmonic I = 1, m = 0. We 
circumvent the use of invariant coordinates/KVFs and instead use the surface averages /i n l of the 
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scalar 2-curvature 2 1Z and Im\l/2 to obtain J\ and higher multipole moments 



A*»(-) :=<«•>- •)"> , (•)■■= j £.dA, (1.3) 

which are well defined, even if the axisymmetry is perturbed and that allow us to benefit from exact 
numerical integration in order to reduce grid size and numerical error significantly. The invariant 
surface integrals fi n ( 2 lZ) , /x n (Im^ 2 ) are related to the horizon spin, mass and higher multipole 
moments by algebraic systems of equations. In principal, the /i n allow to generalize the horizon 
multipole moments through solutions of these systems in the absence of axisymmetry. 

In order to minimize the numerical error of n n ( 2 lZ) , fi n (lm^2) accurate numerical computa- 
tions of the curvature components 2 1Z, lm^2 and the surface triad 2 on the horizon are required. 
The horizon is usually given by h(9, 0) = y/5ijX i X^>, the Cartesian shape function, where X j 
are the Cartesian coordinates at the 2-surface centered at a point inside. Instead of finite differ- 
encing we expand the shape function in terms of a tensor basis to determine Cartesian derivatives 



off the surface, as commonly used in horizon finding algorithms H44I1 . But opposed to Q44] , we 
use another basis, which is easier to implement, and exact numerical integration to determine the 
multipole coefficients of h(9, 4>), where iHl use minimization. 

We apply the new method (in comparison with others) to the dynamical AH of a non- 
axisymmetric BH 3 ringing down to Kerr in a 3+1 simulation, where we follow the evolution of 
spin and mass multipoles until their final Kerr values are reached. 

This paper is organized in the following way. In section fll] we briefly explain the numerical 
methods we use to compute KVFs and aKVFs on AHs. In section [TIT] we deduce formulas from 
the Kerr metric to determine Kerr spin and mass from the area and the 'equatorial' circumference 
or the extrema of the scalar 2-curvature on the horizon and give a new formula which requires the 
surface average ^CTVj and that we also apply to our simulations. In section [IV] we show how 
to use the whole set of [i n to extract the multipole spectrum of an axisymmetric isolated horizon. 
In section |V] we show how to compute the curvature components 2 1Z, and the surface triad 
accurately. In section [VI] we explain the setup and initial data of our 3+1 simulation. During 
the evolution we follow spin, mass and higher multipole moments, compare different methods to 
measure the spin and test their convergence. Notation: Indices i,j, k indicate 3D Cartesian com- 
ponents, indices a, b, c label 2D components on the local horizon grid, letters /, m label spherical 
harmonics. We indicate dimensionless quantities (mass dimension) with a hat, e.g. a = a/m, 
2 TZ = 2 n- A/(8ir), Im* 2 = Im^ 2 • A/ {Air). 



H. SOLVING THE 2D KILLING EQUATION NUMERICALLY 

The IH multipole moments are defined in an invariant coordinate system [|4lll which requires 
the knowledge of the axial KVF on the horizon. Our approach does not explicitly require the 
KVF to extract the IH multipole moments and circumvents the invariant coordinates by using 
the surface averages fi n ( 2 1Z), /i„(Im^ 2 ) which can be easily computed in any coordinate system. 



In statistics fi n is called the nth central moment of the probability distribution of a random variable. 

2 Note that a 'coordinate-induced' surface triad on 'large' coordinate spheres (as for wave extraction via ^4) can 
be easily computed analytically. On the other hand, the coordinate representation of the horizon is a deformed 
2-sphere and the computation of derivatives delicate. 

3 We evolve two puncture initial data with an initially non-axisymmetric common horizon. 
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Nevertheless, in the numerical simulation of section[VI]we want to compare our method and hence 
require the KVF. Therefore, we will briefly explain the techniques we use to solve/approximate 
the Killing equation. 

The induced 2-metric q ab of an spheroid S embedded into Euclidean space admits one rotational 
Killing vector field $ a which is a solution of the Killing equation 

Uq ab = 2 2 D (a ^ b) =0, (2.1) 

where 2 D is the induced covariant derivative on S. The vector field $ a is unique up to a constant. 
For Kerr $ a = d^, where is the Boyer-Lindquist coordinate, this constant is fixed such that 
integral curves have affine length of 2ir, thus G [0; 2%]. 



A. Killing Transport Method 



In order to solve the Killing equation we apply the Killing Transport method Q37D , appendix of 



Il45h . which is explain in this subsection. 

The method can be roughly divided into three steps: 1. determine a single vector of the KVF at 
a point on an arbitrary loop on S, 2. spread this vector throughout the whole surface, 3. normalize 
the whole KVF by normalizing an arbitrary integral curve to have affine length of 2n. The first 
two steps require the Killing transport equation 

c a2 D a <$> b = c a L 2 e ab (2.2) 
c a2 D a (L 2 e bc ) = c a2 R d cba $ d , 

where 2 e ab denotes the Levi-Cevita tensor and 2 R d cba the 2D Riemann tensor. The first equation 
holds, since 2 Dr a & b \ = if is a KVF and since any two-form on S can be expressed as £% ab , 



where L is a function. The second equation follows from the first, see 114511 for details. Therefore 
(12.21) hold for a KVF $ a and the corresponding function L for any vector field c a . 

On the other hand, assume that and L were unknown, pick a loop, e.g. the equator c c , 
[Q = iv/2, 0) of a spherical coordinate system, pick a point, e.g. P, [Q = n/2, = 0) and identify 
c a := d$ 4 , then (12.21) becomes an ODE for the unknown ($i(0), $ 2 (0), £(0)) along c c . This 
defines a linear operator for 3-vectors at P. If we pick three arbitrary, linear independent initial 
vectors at P, transport (12.21 ) them along the loop to P, we obtain a 3 x 3 matrix presentation of this 
operator. Two components of its eigenvector are the KVF at P (1. step), the third is the auxiliary 
function L at P. At next this 3-vector is transported with (12.21 ) along coordinate lines all over 
S, setting c a = or c a = dg respectively (2. step). Where the transportation equation (12.21 ) by 
construction 'conserves' the Killing property. The last step is to normalize the KVF (3. step), 
where we have to solve the ODE d t 9 = 0), <9 4 = $ 2 (6 I , 0) , $q, where the initial vector $g 
is arbitrary, to obtain an integral curve and normalize such that the curve parameter t E [0; 2ir]. 



B. Approximate Killing Vector Fields 



If the spheroid S is slightly deformed, similar to the initial non-axisymmetric AH in our simu- 
lation, no exact solution of (12.11) exists. But one could try to find a 'best match' which minimizes a 



4 The resulting KVF is independent of the initial loop, initial point and curve parameter. 
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certain norm of the l.h.s. of (12.11) on S. Such vector fields are often denoted as approximate Kil ling 
vector fields (aKVF). Opposed to KVFs there is no unique definition of aKVFs. Dreyer et al. 0711 
could show that the Killing transport method is still applicable to yield a 'well matching' aKVF. 
But one has to be aware that the final vector field will not be anymore independent of the particu- 
lar loops of transportation. Although this effect may be negligible for practical applications, e.g. 
[0,0], if the departure from axisymmetry is 'small'. The method has also been used to determine 
aKVFs in binary black hole initial data, see Caudill et al. H39TI . 

We found it useful to adapt the coordinate system on the horizon before applying the Killing 
transport method such that the azimuthal transport revolves the minima of the scalar 2-curvature 5 , 
see appendix El Another approach to find an approximate Killing vector field has been given by 
ll38ll . They use a variational principle to minimize the 'non-symmetric' features of the vector field. 
A similar method can be found in the appendix of [14611 . for an application to a BBH simulation 



see iflfjl . Recently Beetle 
older proposal by Matzner 



[£Zp pointed out that Cook's [13811 approach is closely related to an 
113 611 . where the aKVF is the solution of an eigenvalue problem. An 
outstanding question is still the normalization of these aKVFs. An interesting new idea has been 
given in the appendix of |46J], where the aKVF is normalized to a particular surface integral instead 
of a single integral curve. 

In our approach these difficulties do not appear because no KVF/aKVF is explicitly required to 
represent the axisymmetry/perturbed axisymmetry. Instead we compute the invariant surface av- 
erages fj, n which exist in any case and from those compute the IH multipole moments/ generalized 
IH multipole moments through the algebraic system linking the two sets of invariants, subsection 

nYAi 



C. Coordinate Vector Fields 

If the coordinates are conveniently adapted to the metric manifold, the coordinate vectors can 
automatically generate symmetries (if existing), such as the Boyer-Lindquist coordinate vectors d t 
and in a Kerr spacetime. This is also the case for the adapted spherical coordinates (# aS c, 0asc)> 
see appendix[Bl and the particular initial setup we chose in our simulations 6 . Then the coordinate 
vector field 

*Lc = 9^ sc , (2.3) 

is a good approximation to the KVF and we can estimate the spin J[$ J ] ps ^/[$i sc ] with (11.11) . see 
the application in section IVIl 

Similarly [3 H^] use the three rotational Killing vectors of Euclidean space in Cartesian coor- 



dinates 

$W = {x k -C k )& k ,j = l,2,3, (2.4) 
where e l \S pk = e ljfc is the flat space Levi-Cevita tensor and a point inside S, to define a Eu- 
clidean spin vector (J[$£ ] ], J[$x], J[#x]) and together with £□} to estimate J[¥\ w J[¥ cc ], 
where J[$cc] denotes the Euclidean norm of this vector which allows them to study the spin pre- 
cession in a BBH inspiral and to estimate the final spin after merger. Referring to 0] this Euclidean 



5 An spheroid has two minima of the scalar 2-curvature which coincide with the minima of the KVF, given by the 
symmetry axis of the body. 

6 In general this is not the case and the correct solution of the Killing equation has to be found. In our case coordinate 
vector fields are very useful for the comparison of section [VT1 
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spin vector reproduces the Bowen- York spin parameters of the conformally fiat initial data and for 
the final black hole | J[&] — ^ 1 as in our simulations. 



III. INVARIANTS OF THE HORIZON IN KERR 

Before we go into the details of how the surface averages fx n are linked to the IH multipole 
moments in the next section [IV] we want to remind that the mass M z Kerr and angular momentum 
jKerr mu itip i e moments of Kerr M z Kerr + iJ / Kerr = m(iJ/m) 1 are uniquely given by Kerr spin J 
and mass m. In this section we will review the analytic formulas necessary to extract Kerr spin 
and mass from an AH and give a new formula which we apply in our simulations. 

In many numerical simulations Kerr spin and mass (J, m) are being computed from the 'equa- 
torial' circumference 7 and the area (L(c e ), A) of the BH surface, see [40]. A more recent ap- 



proach is to use an extremum of the scalar 2-curvature and the area ( 2 7£ ext , A), see ll46ll . 
Each of these invariant pairs uniquely determines a Kerr spacetime and is related to the other 
through the Kerr metric such that we are free to choose the numerically most convenient one. 
In order to benefit from exact numerical integration and to avoid interpolation on the hori- 
zon we chose the invariants (yU 2 ( 2 7?.), A), see (11.31) . The explicit algebraic expressions relating 
J <-> L(c c ) <-» 2 1Z CJ it <-* yLt 2 ( 2 7?.) (<-> /i 2 (Im^ 2 )) are derived in the following. 
Any axisymmetric 2-metric q a b can be put in the compact form 

d <? = 4- ( -7^ d X 2 + fixW 2 ). (3.1) 



47T V/(X) 



For the 2-surface of a Kerr black hole f(x), see [48], is given by 



/(*) = - \ 2( f 2V X--=cose, (3.2) 
1 - (3 2 (1 - x 2 ) 

where (3 G [0; 1/V2] is called the Kerr distortion parameter and (8, (f>) are the Boyer-Lindquist 
spherical coordinates. The distortion parameter (3 is related to the more familiar dimensionless 
spin parameter a = a/m = J/m 2 by 

P 2 = I ( 1 - Vl^ 2 ] = ^-77, (3.3) 



to Kerr spin J and mass m 8 by 




2 V J c 



A (3 




*Jl-fr 871 2 Y 4tt(1 - /3 2 ) 



(3.4) 



Smarr H48H pointed out the analog of the surface of rotating material bodies to the black hole 
horizon, where the equatorial circumference increases as the body spins up. The equatorial cir- 
cumference for the Kerr horizon is given by integrating (13.11) along the maximum of 2 1Z which is 



7 This is the curve c c along the maximum of 2 1Z in Kerr. 

8 For completeness note that mj rr = i? a reai/2 is the irreducible mass and i? arC ai = y/A/ (4tt) the areal radius. 
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the curve (x = 0, 4>), 



f / A / Att 

L{Cc) = t V4^ /(x = 0) ^ = VTT^ = 47rm " (3 ' 5) 

For the numerical application in arbitrary coordinates this is practical, if the curve c c is known to 
overlap with a coordinate line. If this is not the case the extrema of 2 1Z are an appealing alternative, 
see tlalfig]- The scalar 2-curvature of q a b (13.11) is 

2 K = -~f"(x) - 2 n = - 1 -f(x), 0.6) 

with extrema at Xmin = 1; — 1 , Xmax = 0. We obtain 

^max ~ ^ ) 1 4/3 . (3.7) 

(l - (3 2 y 



A. An invariant surface integral in Kerr 

If the scalar 2-curvature (or alternatively Re\l/ 2 , since Re\I/ 2 = — 1 /4 2 72. for Kerr) has been 
computed on a finite grid, interpolation is required to obtain the extrema. This is not necessary if 
the following surface integrals are employed 

^n):=^n)- 2 n) 2 y ( 2 n) : = jf^dA. (3.8) 

Moreover, the numerical error of /i 2 ( 2 7?.) benefits from averaging over all points on the grid and 
exact numerical integration can be used. With the normalization of (13.61 ) the average < 2 7Z > gr id= 
1 + fnum, where e num is the numerical error, for any 2-metric computed on a finite grid on S 
according to the Gauss-Bonnet theorem. For Kerr the integral appearing in (13.81 ) is taken over a 
rational function in x- We obtain 

2 ~ _ -15 - 70c 2 + 128c 4 + 70c 6 + 15c 8 3(1 + c 2 ) 4 arctan(c) 

/M k>) — 1 — , (3.9) 

80(1 + c z ) 16 c 

where c is defined in (13.41) . In our simulations we compute the surface average /i 2 ( 2 7£.) numer- 
ically and solve (13.91) for the Kerr c. Kerr spin and mass are then given by J = A/(8ir)c and 
m 2 = A(l + c 2 )/(167r) (I3.4I ). For the numerical application in section Ivn the Kerr spin devi- 
ates significantly from the IH spin during the initial phase but the 'non-Kerr' features are radiated 
during the evolution and finally vanish below the numerical error. 

Note that we could similarly use any fi n ( 2 7Z), n > 2 or ^ n (Im^ 2 ), n > 1 to compute c for 
Kerr. In that case Im^ 2 = — \g"{x)i dix) '■= c(i+d 2 x 2 ) ' see anc ^ some algebra. It follows that 
^!(Im^ 2 ) = and /ii(x • Im^ 2 ) = c. The explicit appearance of the Boyer-Lindquist coordinate 
(X = cos 9) is inconvenient for the numerical application. For /x 2 (Im\^ 2 ) we obtain an expression 
similar to dH which is // 2 (M 2 ) = ^+i7oe^g^p!+i5eg + 3(i+p* arctan^ Tq extract more 

information than the Kerr c we have to consider the whole set of fj, n and follow the procedure 
explained in the next section. 
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IV. INVARIANTS OF AXISYMMETRIC ISOLATED HORIZONS 



For the calculations in the last section to be reasonable when applied to an AH found in a 
numerical simulation, we had to assume that the detected 2-surface was in a slice of Kerr. We 
relax this condition and allow the spacetime to be dynamical in the vicinity of the horizon which 



we assume to be an axisymmetric isolated horizon (IH) 11321 . 13311 . On the horizon in Kerr all 
multipole moments are necessarily given by spin and mass, therefore higher moments contain no 
extra information. This is in general not the case on an axisymmetric IH, where an infinite set of 



independent multipole moments permits more complexity, see 114 111 . 



Ashtekar et al. 114 111 exploit the axisymmetry to define an invariant coordinate system 0) for 
which the 2-metric has the form (13.11) . is the KVF and the (zonal) harmonics {Y l0 (x)} represent 
an orfhonormal basis <f> s Y 10 (x)Y l '° (x)dA = which they use to define the dimensionless IH 

mass 1 1 and angular momentum L\ multipole moments 

/,:= <f l/4 2 n( X )Y m (x)dA, L r .= - (flm^ 2 ( X )Y l0 (x)dA. (4.1) 

Js J s 

. oo . oo 

2 U{x) = 4--J Y l0 ( X ) , M 2 ( X ) = ~ Yl ^ Yl °M ■ ^ 

l=0 l=0 

On IHs without matter fields (like in Kerr) the Weyl scalar ty 2 is invariant and Re^ 2 = — 1/4 2 72.. 

Note that for Kerr J ■ 8n/A = c = ^1/(3tt) Li and for an IH J[&] -8tt/A= >/1/(3tt) Li, 
where & is the KVF corresponding to (x, 0) and J[&] given by (11.11) . Therefore, the curva- 
ture component lm^2 is sometimes called rotational Weyl scalar and the Li angular momentum 
multipole moments, all vanish in the absence of spin. 

The invariants I h Li are subject to certain algebraic constraints such that I = y/n (Gauss- 
Bonnet), that the mass dipole 1\ and the angular momentum monopole L vanish 9 . If the 2-metric 
(13.11) admits a reflection symmetry as for Kerr f(x) — f(— x)> see (13.21) . all odd and even Li 
vanish, too. 



A. The invariants \i n on axisymmetric isolated horizons 

In analogy to the method explained in subsection IIII Al for Kerr, where we gave the formula 
(13.81) to compute the Kerr c from the surface average fi 2 { 2 7Z), we would like to relate the invariants 
fj, n ( 2 7Z), /i n (Im^ 2 ) (|1-3I) . which are numerically easy to obtain in any coordinate system, to the IH 
multipole moments (14.11) which would require the invariant coordinates for a direct computation 
of the integrals (14.11) (as for example being done in ll50ll ). 

We obtain the algebraic relations between the ^ n ( 2 TZ), /i n (Im^ 2 ) and the I n , L n by inserting 



9 Therefore, the invariant coordinates are sometimes called 'center of mass frame' of the IH. 
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(S3 into (TO) 



( 2 7t) = /(l-2gj,F w ( X ) j\ , n = 2,3,...,n max , (4.3) 




nM 2 ) = ( \0 + YLY m (x)\ ), n = 2,3,...,n* (4-4) 



where we assume that 2 TZ, Im^ 2 are given by finite sets of multipole moments up to Z^ ax , I 
We obtain 



L 

max' "max' 

10 



n (lm$ 2 ) = E (j/ )(^-)^ max <( r "°)^ max > ' ™ = 2,3,...,n max , (4.6) 



l A 'imaxl =rl 

where -K"z max = fc 2 , &z max ) is a multi-index of length / max , ( fc * ) is the multinomial coef- 
ficient and (IJ*w ((y-Q)^w) = (f 1 )*i(J 2 )fe...((y 10 )fci(y 20 ) fc a„.) . The integers n max ,n max 
match the numbers of non-trivial I n , L n given by the algebraic constraints mentioned earlier and 
^max i ^max- The coefficients ((Y-°) Kl ™** ) are integrals over products of (zonal) spherical harmon- 
ics. They are given by the associated Clebsch-Gordan coefficients and higher order generaliza- 
tions. 

Consider the following example. In a simulation of a perturbed Kerr spacetime we locate the 
AH and compute the surface integrals Li n ( 2 7V) ,/i n (Im^ 2 ) (11-31) numerically to n max = 6 n to 
equate them with the r.h.s. of (14.31) . where we assume that the 2-surface is a cross-section of an IH 
with reflection- and axisymmetric 2-metric. Then the algebraic systems (14.31) , (14.41 ) become 

// n ( 2 ft) = / (l - 2(v^r 00 + Yl IiY m ) + o)j Yn = 2,3,4,5,6 (4.7) 

\ V £=2,4,6,8 / / 

/i„(m 2 ) = ^J2LiY l0 + O L ^ ^,n = 2,4,6, (4.8) 

which we solve for I 2 , I4, Iq, Ig, Oj and L±, L 3 , Ol, where Oj, Ol are constants accounting 
for the truncation of the expansions. Since we simulate a perturbed Kerr spacetime, we pick the 
solution that is real and for which — 7 2 > I 4 > —I 6 > Oj and L x > —L 3 > Ol holds as for Kerr. 

In analogy to electro dynamics dimensionfull factors can be added to attribute a physical inter- 
pretation to the Ii, Li, see JIj]]. To obtain the spin we need 



Here the indices I, Lin Z^ ax , Z^ ax , n^ ax , n£ lax are omitted. 

Formally the solutions of the algebraic systems depend on n max . It determines the number of multipole moments 
we can resolve l max and is limited by the numerical noise. In pratice the solution for lower Z max does not change 
as we go to higher n max . 
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j > = Vt^" <4 - 9) 

The equation J[&] = J\ holds if $ J is the KVF corresponding to the invariant coordinates (%, 0). 

The surface integrals /i n are well defined even in the absence of axisymmetric and allow to 
extend the concept of IH multipole moments by adding the m/fl harmonics in the expansions of 
2 1Z and Im^ 2 on the r.h.s. of 14.91 Nevertheless, for the evolution of the non- axisymmetric initial 
data studied in|VT]we assume that the contribution of ode/even mass/angular momentum multipole 
moments (reflection symmetry), higher harmonics as well asm ^ harmonics is small and can 
be accounted for through Oi, Ol- We do not further investigate the possibility of generalized 
multipole moments. Our approach aims at numerical convenience and is flexible enough to extract, 



in principle, other invariants like the generalized multipole moments proposed by Owen [51] who 
considers the eigenfunctions of the intrinsic Laplacian on the horizon. 

V. ACCURATE COMPUTATION OF 2 TZ, f 2 ON THE AH 

In this section we will show how to compute the curvature components 2 1Z and ^ 2 accurately, 
where we assume that the 3+1 evolution variables 12 extrinsic 3-curvature K^, 3-metric 7^ (to- 
gether with diKj k , di'jjk, didj^kk') and the horizon coordinate shape X'-* are given on a Cartesian 
grid. The accurate calculation of curvature components on a deformed 2-sphere in a Cauchy 
slice is a common problem in numerical relativity which appears in horizon finding algorithms. 
Various methods have been tried to discretize the necessary spatial derivatives djh, didjh by finite 
differencing, finite element, pseudo-spectral and spectral methods, using squared (9, 6) grids or 



multipatch grids, for a review see [520. Our approach is motivated by the work of H44H . There a 
spectral decomposition of the coordinate shape function h(9, 0) is being used to compute Cartesian 
derivatives. The 1st derivatives djh are necessary to obtain a surface triad {s\ u\ v k } (required 
to compute the Weyl scalars) and the 2nd derivatives didjh to obtain the extrinsic 2-curvature 2 Kij 
of S embedded into the Cauchy slice (additionally required to compute the scalar 2-curvature). 

If we parametrize the AH with spherical coordinates, the embedding (9, 0) into the Cartesian 
grid is 

X 3 {6,<f)) = h{6,<j))n 3 +C j , (5.1) 

where is a coordinate location inside the horizon (for example the coordinate centroid), n J the 
Cartesian radial unit vector n- 7 = -x\ r = yj b^xi and x- 7 are Cartesian coordinates. 

A. Spectral decomposition 

To compute spatial derivatives one could decompose h(9, 0) into 

h{d^) = Y,Y,^ lmYlm ( e ^^ (5 - 2) 

Z=0 m=l 



where [h] are the expansion coefficients and Y the spherical harmonics. The evaluation of 



djY lm (9, 0) would require the Jacobian to transform between spherical and Cartesian coordinates. 



12 They can be easily assembled from the BSSN evolution variables. 
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This is inconvenient for the numerical application, since the Jacobian is singular at the spherical 
coordinate poles. 



Therefore, Il44h take a tensor basis which is build of the radial unit vector n l (x^) = x J /r and 
thus defined in Cartesian coordinates (and easily parametrized with any other local coordinate 
system on the 2-surface, e.g. spherical n l (9, <fi) = (sin 9 cos 0, sin 9 sin 0, cos 9) or stereographic 
coordinates n j (u, v) = (2m, 2v, u 2 + v 2 — 1)/(1 + u 2 + v 2 )), 

Irnax 

h = J2[h] Kl N Kl , (5.3) 

1=0 

where K[ is again a multi-index of length I, N Kl = n kl nk 2 ■ ■ ■ n kl is the vector product of unit 
vectors and the location-independent coefficients [h] Kl are symmetric tracefree tensors (STF), the 



notation is adapted from [153( 1 . If the STFs are known, they can be translated to obtain the expansion 



(15.31) . for how to [h] lm <-> [h] Kl see [14411 . The partial derivative of the tensor product djNx t consists 
of the derivatives dirij = (5ij — nirij)/r. In detail the implementation of the STF tensors and its 
derivatives is a bit cumber somely but straight forward. 

We use another basis of the harmonics instead {5ijn l N^) 1 , where A/"- 7 is a constant complex 



Euclidean null vector (MjJ\f j ) = 0, N j ^ 0, see Sec. 11.5.1., Vol.11 D54[] or [55J]. The expres- 
sion (njM^) 1 is a homogeneous harmonic polynomial of Euclidean space of order I, therefore 
Afl a t (nj.APy — 0. The radial vector -nP defines a restriction of the polynomial to the unit sphere 
x l xH{j = 1. It is known that such restrictions are eigenfunctions of the Laplacian of the induced 
metric (this applies to any embedding of S 2 into Euclidean space, e.g. an ellipsoid). On the unit 
sphere this implies AofajNi) 1 = 1(1 + l)(rijj\f^) 1 , where A D is the Laplacian of the standard 
spherical 2-metric. This holds for any null vector J\fi. In order to span each /-eigenspace of A D 
with 2/ + 1 linear independent eigenfunctions we define a list of null vectors 

2ix 

■^[im] = (ism(mai),icos(mai), 1) , a t = ^ ^ m = -/,•••,/ , (5.4) 

where the roots of unity have been used such that the A/j], , have the Euclidean norm A//.AP = 
_| e *2i+r | 2 _|_ i Now we can define the new basis $' m := (njAf^) 1 and decompose h into 

'max I 

^EE^W- (5 - 5) 

1=0 m=l 

The <3>' m , m = —I, ■ ■ ■ ,1 are not orthogonal in each /-eigenspace but across different eigenspaces. 
They are related to the standard basis by 

i 

m'=—l 
i ' ~ylm' 

(him _ \ ^ im'mai sr 7 \ 

21 + 1 2^ B lm> e ' V''> 

m'=-l 



R im _ ( nm l (l + m)\(l-m)\ 
B ~ UV 4tt(2Z + 1) 
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and we can transform the coefficients \h] lm <-> [h] l ff. Derivatives of the new basis 13 are given by 

d k $ lm = ( nj Af s ) l -H {d knj M j ) (5.8) 
d k & m = (nj^y-H ~ (A4 - n knj M j ) , (5.9) 



and similarly for higher derivatives QjC>,$ 



B. Surface triad 



lm 



Now we have the Cartesian derivatives djh, djdih at hand and are able to compute the outward 
pointing surface normal s J = ^ k s k 

8 . = - djh), A = l/yj^^m-dih^rij-djh). (5.10) 

In order to complete the surface triad {s\ u\ v k } we set w- 7 = , 1 ^=dpX^ and v k = 

y/likdeX i d e X k 

e^ k SiUj, where e l ^ k = ||7||~ 1 / 2 [123] l - ?fe is the spatial Levi-Civita tensor and [123] yfc the pure alter- 
nating symbol. 



C. Extrinsic and intrinsic 2-Curvature 

The extrinsic 2-curvature 2 Ky of S embedded into the Cauchy slice is given by 

2 K tj = D iSj - s lS k D kSj , (5.11) 

where the second derivatives djd k h are required and the Christoffel symbols associated with the 
3-metric to compute the 3-covariant derivative Dj. Then the intrinsic 2-curvature 2 1Z is given by 
Gauss' theorema egregium 

2 n = TZ- 2R ij s { s j + X 2 - 2 K lj 2 K i3 , (5.12) 

where 2 /C = 2 Kijq % '-> and g y = ^ — sV is the induced 2-metric in Cartesian components (also 
required to raise the indexes of 2 Kij in the last summand on the r.h.s. of (15.121) ) and Rij, TZ are the 
3-dimensional Ricci tensor and scalar. 

D. Area Element 

The computation of surface integrals on the AH requires the area element dA = y/det q ab dddcj), 
where we need the induced 2-metric in local coordinates 

q ab = d a X^d b X k ljk , (5.13) 



here X J has been defined in (15.11) . for an alternative see appendix of 114411 



13 Here we omit the subscripts A/^ m j — ► A/" 3 for simplicity 
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E. ^ 2 and other Weyl scalars 



To obtain mass and angular momentum multipoles (14.11) an accurate computation of ^2, given 
the 3+1 evolution variables, is requried. Additionally, we want to follow the constraints \[> = 
and \&i = which hold for Kerr and on IHs ll3~2ll in the simulation of section [VO The electric 
and magnetic parts of the Weyl tensor Cij k i w.r.t. time-like normal of the Cauchy slice are 



Eij = 
B i:j = 



-Cijkin k n l = —Rij + Ki k K k j 



KK, 



1 j ■ 



-Ei kl D k K u • 



(5.14) 
(5.15) 



— * Cijkifi h' = —ti uk^ij 
We further project E^, Bij onto the surface triad {s\u j , v k } and obtain the Weyl scalars, see 



*2 

*o 



--(E jk -iB jk )s j s k , 
-(E jk - iB jk )m>m k , 
-—={E jk - iB jk )m J s k , 



(5.16) 
(5.17) 
(5.18) 



where m? = -^(u j — iv j ). 

We monitor the dynamics of the AH during the evolution in section [VI] by computing the 
dimensionless surface integrals 



4>o 



|*i I dA, ^ 2 



s 



— d> 4Re * 2 dA + 1 
8tt ^ 



(5.19) 



which vanish for a MOTS in a slice of Kerr or an IH. 



VI. NUMERICAL EVOLUTION AND INITIAL DATA 



In order to test and compare the new techniques we applied them to the dynamical AH of a 
non-axisymmetric spinning BH in a 3+1 simulation ringing down to Kerr which as been carried 
out using the CCATIE code 111 in . This is a 3D finite differencing code based on the Cactus Com- 
putational Toolkit 115 811 . The CCATIE code provides a collection of modules (thorns) which allow 
us to use puncture initial dat a Il59ll w ith the TwoPunctures thorn 11601 . to do time evolution using 
the BSSN evolution system 116 ll . l62l 16311 . to set proper gauge conditions (where we used 1+log 
slicing and a hyperbolic gamma-driver condition stemming from Il64ll but with advection terms 
llllll ). to successively refine the Cartesian mesh with several nested static boxes around the AH 
(where we used the Carpet AMR driver I165ll ) and to locate the horizon every few time steps during 
the evolution ll66ll . The horizon finding thorn provides the shape function h(6, <p) which is being 
used by a separate thorn to interpolate (4th-order Lagrange) all necessary 3+1 evolution variables 
onto the spherical grid, to accurately compute the curvature components 2 1Z, Im$2 at the horizon 
(see section [V]) and, finally, to determine the associated quasi-local IH multipole moments using 
the surface integrals jjL n (11.31) . 



A. Initial Data and Grid Parameters 

In order to model the common horizon after the coalescence of an arbitrarily aligned BBH 
system we chose as a non-trivial initial configuration a misaligned spinning puncture with a nearby 
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smaller non-spinning companion puncture, where the common horizon is already present on the 
initial slice. The Bowen-York parameters of the first puncture are mi = 0.8M, |si| = 0.3M 2 
with orientation (9 S1 = 0.6, </> Sl = 0.4) in the Cartesian grid. And for the second puncture we set 

m 2 = 0.2M, s 2 = 0. 

The evolution is being carried out using the method of lines with 4th-order Runge-Kutta time 
integrator and 4th-order centered stencils for spatial differentiation with the Cartesian grid res- 
olutions Ax = 0.048M, 0.035M, 0.025, 0.02M (finest AMR resolutions). To determine the 
KVF/aKVF we use the Killing transport method III Al with 2nd-order centered stencils for dif- 
ferentiation and a 2nd-order Runge-Kutta integrator, see [37] for more details. To compute spatial 
derivatives of the shape function h(9, <p) we use its decomposition into Spherical Harmonics where 
the spectral resolution is fixed to Z max = 10. To compute the surface averages fj, n we use an exact 
integration scheme, see appendix lAl and fix n max = 6. For every Cartesian resolution we use 
three different spherical horizon grid resolutions N e x = N s = 480, 1104, 4900, where N$ is 
the total number of grid points on the surface and = 2(N e + 1). The horizon finder is using 
a projective 6-patch grid Il66ll with approximately the same number of points as on the spherical 
grid. 



B. Numerical Evolution 

1. Monitoring the Isolation Constraints 

To monitor the dynamics on the horizon we computed the surface integrals (15.191 ) shown in 
figure[T](for Kerr ^0,1,2 = 0). On the left we see the typical exponentially damped oscillation of the 
radiative Weyl scalars \l/ , which are (after an initial burst ^ ,i ^ 1) given by a superposition 
of several quasinormal-modes, predominately / = 2 modes, that have been excited by the specific 
initial data. As a fit to the ring-down profile of ^ we obtain the frequency u)& t « 0.355 + O.O882, 



in agreement with the / = 2-mode frequencies c^ =2 . mri , see H67H . which are cj 2 -2o ~ 0.34 + 
0.089i, U220 « 0.36 + 0.089i, • • • for the case J = 0.3, m = 1.035. After around t > 90M 
the perturbations are to weak to be further resolved limited by the total numerical error, which 
we downsize by increasing the Cartesian grid resolution, see figure [Hon the right, in order to see 
the dynamics below ^ < 10~ 5 . For Ax = 0.035 (black and orange) we computed xjj for two 
different spherical resolutions to show that the total error of ^0 (and similar for surface integrals of 
other curvature components) is almost independent of the spherical resolution due to the spectral 
methods involved. 



2. Evolution and Convergence of the Invariants fi n 

In figure[2]we see the exponentially damped oscillation of the fi n as they ring-down to their final 
Kerr value. On the right it is shown how the time averages of ^i^H) (120M-200M, straight black 
lines) converge with the expected 4th-order (4.01) as the Cartesian grid resolution increases after 
the oscillations have settled down. Apparently, the error of /i 2 ( 2 7?.) does not converge uniformly 
but the effect flattens out as the Cartesian resolution increases. 
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FIG. 1: Left: time evolution of dimensionless surface averaged Weyl scalars "00,1,2, Right: time evolution 
of V>o for 3 different Cartesian resolutions 
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FIG. 2: Left: time evolution of the surface averages /i2,3,4,5( 2 ^) over powers of the scalar curvature of the 
horizon, Right: time evolution of /i2( 2 ^) for 3 different Cartesian resolutions and time averages (straight 
lines) between 120M — 200M for each resolution 



3. Evolution of Mass and Angular Momentum Multipole Moments 

From the \i n we compute the IH multipole moments J/, Li corresponding to an reflection and 
axisymmetric horizon by solving the algebraic system (14.71) . where Oj, Ol account for all higher, 
non-axi symmetric and non-reflection symmetric multipole moments. It is apparent in figure[3]that 
these multipole moments are quickly radiated t < 30M, leaving the horizon almost reflection and 
axisymmetric but still oscillating. Interestingly, the dimensionless IH spin L\ is almost constant 
during the evolution, as the horizon area (not plotted, A « An ■ 2.05 2 M 2 ). 

4. Spin Evolution and Comparison with other methods 

In figure |4] we see the comparison between the various spin measures and their convergence. 
We have 

1. Ji — v4/Vl927r 3 Li (red) computed from the fi n , (14.71) . assuming an axisymmetric IH, 
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FIG. 3: Left: time evolution of mass 1\, Right: and angular momentum multipole moments Li given as a 
solution of the algebraic system ( 14.71 ) for the fi n up to n max = 6; Oj, Ol account for all higher multipole 
moments 
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FIG. 4: Top Left: time evolution of spins given by the Killing transport aKVF <J>k t , the coordinate vector 
fields $aso ^cc and the Ken - spin computed from ^(^IZ); Top Right: Zoom of 'Top Left' together with 
angular momentum dipole J\ = A/V 1927r 3 Li (red) computed from fi2, fJ>4, of Im^; Bottom Left: 
convergence of J[<3?kt] varying number of spherical grid points Ns; Bottom Right: convergence of J[$ asc ], 
J(/i2( 2 7£), A) varying Cartesian resolution Ax 
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2. J[$kt] (blue) computed from the Killing transport KVF/aKVF $ kt , (f2T2b . (fTTTT) . assuming 
an axisymmetric IH, 

3. Jf^cc] (light green), J[$ aS c] (dark green) given by the coordinate vector fields $ cc (Cartesian 
coordinates), (12.41) . $ asc (adapted spherical coordinates), (12.31) . assuming 'small' coordinate 
distortions, 

4. J = J(/i 2 ( 2 7?.), A) (brown) computed from yU 2 ( 2 7£), (13.81 ), assuming a Kerr horizon. 

After a short initial bust all methods yield nearly the same spin value, which stays constant dur- 
ing the evolution; except J(yU 2 ( 2 7£), A) (brown) which oscillates with the quasinormal frequency. 
During this phase the horizon seems to be best modelled assuming an axisymmetric dynamical 
horizon but not Kerr. We chose the numerical setup such that the coordinate distortions are small 
and J[$as C ], «/[$cc] overlap with the invariant measure J\. This is in general not the case in a full 
BBH simulation and these methods should be used with care. 

In figure H] (bottom right) we see the expected 4th-order convergence (w.r.t. Cartesian grid) of 
^[•^asc], ^[$cc] an d J(/i 2 ( 2 7^), ^4) towards 0.3M 2 . The convergence of J\ is not shown explicitly. 
It is a smooth function of the fi n (convergence shown above) and converges therefore at the same 
rate. On the other hand J[<i\t] converges at 2nd-order (w.r.t. the spherical grid), figure |4] (bottom 
left) 14 , because the Killing transport method requires finite differencing on the horizon grid to 
determine $£ t . 

VII. CONCLUSION 

The dominant part of the gravitational radiation at Scri is contained in the quadrupole moment 
of ^ which is in practice extracted at 'large' coordinate spheres around the source in numerical 
simulations. Similarly, the dipole moment of the rotational Weyl scalar Im\[> 2 encodes the quasi- 
local angular momentum measured at the apparent horizon in the presents of axisymmetry. The 
local coordinates on the horizon are in general distorted and a solution of the Killing equation 
is required to determine an invariant coordinates system in which the multipole moments can be 
computed. 

It is involved to determine the Killing vector field, in particular, to find a convenient approxi- 
mant in case the axisymmetry is perturbed. We have shown a new method to extract the horizon 
multipole moments using coordinate invariant surface integrals fi n from which we deduce the mul- 
tipole moments as a solution of an algebraic system. In case of an axisymmetric IH the angular 
momentum dipole J\ is equal to the spin J[$] given by a solution of the Killing equation $ J in 
agreement with our simulations. Interestingly, the spin of the aKVF $ kt (given by the Killing 
transport method) and the angular momentum dipole moment Ji(/i n (Im\E' 2 ) , A) (given by the /i n ) 
agree even in the absence of axisymmetry. 

There seems to be a dynamical phase of the horizon in which it is better modelled by an ax- 
isymmetric dynamical horizon and not with Kerr. Nevertheless, after the horizon is settled the 
Kerr formula is valid. Then the computation of the Kerr spin using the surface average /i 2 ( 2 7?.) 
(or /i 2 (Im\l> 2 )) is sensible and numerically more convenient than using the horizon circumference. 
The deviations from Kerr oscillate in agreement with black hole perturbation theory, until they are 



Note that the low resolution Ns = 480 (light blue) is to coarse to be in the convergence regime. 



17 



no more resolvable due to numerical errors. Then the dipole moment of the rotational Weyl scalar 
agrees with the Kerr spin and the ji n take their final Kerr values. 

We have shown how to use spectral methods, in a 3+1 finite differencing code, to accurately 
compute curvature components at the horizon and to extract spin and other multipole moments 
saving computational costs. These techniques, in particular, the non-standard basis of spherical 
harmonics and the exact integration scheme, should be considered for wave extraction on coordi- 
nate spheres or constant mean curvature spheres |68i 69ll . 
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APPENDIX A: EXACT INTEGRATION SCHEMES FOR SPHERICAL GRIDS 

It is well know that the equation 

/ f(x)w(x)dx = y]wif(xj) , (Al) 
Ja i=i 

holds exactly, where w(x) is called the weight function, if f(x) is a polynomial of degree less than 
2N and the weights W{ and abscissas X{ are chosen in accordance with the orthogonal basis of 
polynomials on [a, b] defined by the scalar product < f\g >:— f(x)q(x)w(x) dx, because there 
are 2N degrees of freedom to make both sides of (IA1I) match, see for example 117 Oil . 

For the integration with w(x) = 1 on the circle a = b, the 'correct' weights and abscissas 
are particularly simple. They are N equi-distant points with equal weights. This can not be 
generalized for the integration on the 2-sphere 

N s 

f(?,y)dA=y2wif(xi,yi), (AT) 

for arbitrary Ns, because the number of uniform grid structures is finite Ns = 4, 6, 8, 12, 20, 
corresponding to the faces of the platonic solids. Since this is a 2D integration, we have 3Ns 
degrees of freedom in the sum on the r.h.s. of (IA2I) and (Z max + l) 2 spherical harmonics of degree 
< Zmax- This means if f(x,y) was given by an expansion up to / max , we needed at least Ns = 
(Imax + l) 2 /3 points to make (IA2I) hold. Lets say f(x) was given by an expansion of (7 + l) 2 — 4 
spherical harmonics, then the integration (IA2I) on an icosahedral grid N s = 20 with equal weights 
would be exact. There is an extensive body of work on the problem of optimal integration schemes 



for N s > 20 (cubature problem), see for example 117 111 . 
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There are less optimal compromises available, which require much more points than (/ max + 
l) 2 /3, but which are defined on regular spherical (9, <p) grids. For example the Gauss- 
Legendre/Gauss scheme, where the integration along each interval [—1, 1], [0; 2tt] is a Gaussian 
quadrature 

f / f(xA)dxd<P = TTw?wff( X i,h), (A3) 
Jo J -i ttU 

where again \ — cos N s — Nq x and = 2Ng. 

As before the 0-integration is a Gaussian quadrature for <frj = 2%(j — l)/N^, j = 1, A^ and 
equal weights Wj = 2-n/N^, the ^-integration (in that case called Gauss-Legendre quadrature) for 
Xi being the roots of the Legendre polynomials (according to the weight function w(x) = 1)- The 
corresponding weights wf can be found in e.g. IT72il . This method is exact for polynomials of 
degree less than 2N e (less than y/2N s < y/3N s ). 



An alternative integration scheme has been found by Q73D 15 . There the integration grid is a 
standard equi-angular (0,0) grid, 9j = (j — 1/2)% /N g (staggered) and the computation of the 
roots of the Legendre polynomials not necessary. The weights for even/odd N e are given by 

Ne/2-l 

Wj = A/N e jsm((2k + l)ej) , iV e even, (A4) 

k=0 

I , (7V 9 -l)/2-l \ 

w * = 4/Ne [ 2N e sin(iYe ' dj) + ^TTT sin m + 1)9 j) J ' Ne odd ' (A5) 

which allows for exact integration of harmonics of order less than N e /2 (less than 
y/lJWs < \ f 2N^ < V3As)- Then equation dA3]) becomes 

p2tt j-* N 6 N <l> 

& j f{9,<f>)sin0ded(t> = ^J2w e i wff{e i ,<j> j )sm6 j . (A6) 
Jo Jo i=1 j=1 

A small summarizing example: for the total of Ns = 512, Ng x — 16 x 32 the cubature limit 
is at 39 ps V3 • 512 = \/3Ns, for the Gauss/Gauss-Legendre scheme we get / max < 32 = 2N e and 



for the scheme of 117311 we have / max < 8 = N e /2 (we get almost the same limit on an icosahedral 



grid 16 with only N 5 = 20, where l max < 8 « V3 ■ 20). 



APPENDIX B: ADAPTED SPHERICAL COORDINATES 



Before solving the 2D Killing equation on a sphere it is useful to have the 2-metric or the hori- 
zon shape in a convenient coordinate representation, which is 'roughly' adapted to the axisymme- 
try. Such that the poles of the spherical coordinates system agree with the two minima of the scalar 
2-curvature. We assume 2 1Z(9, 0) to be given on a spherical coordinate system (9, <fi), where the 



The authors make use of the fact that the points Xj = cos @j (although not the zeros of the Legendre polynomials 
on [1; —1]) are the zeros of the Chebyshev polynomials of the 1st kind. 

Therefore, if one is only interested in the first coefficients of a smooth function on the sphere up to i max = 6, an 
icosahedral grid with equal weights would be a good choice. 
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FIG. 5: parametrization of the unit-sphere with a shifted spherical coordinate system 



two minima are already in the xz -plane symmetric to the x-axis at iV- 7 = (sin 9 min , 0, cos 6> min ) and 
S j = (sin0 min , 0, — cos^min), see figure |5] This can always be accomplished by a simple Euler 
rotation. 

In order to obtain the adapted spherical coordinates system (9', 0'), we have to shift the Carte- 
sian z-axis along the a>axis by the amount d := sin 9 min . This is being done by 

0) = r'{6', 0') n'^9', <f>') + d- (1, 0, 0) , (Bl) 

where rij(9, 0) = (cos sin 9, sin sin 9, cos 9), n'^{9' , 0') = (cos 0' sin sin 0' sin cos 6 1 ') are 
the radial unit vectors in the corresponding coordinate system. The distance r'(9', 0') is given by 

r'(6', 0') = Jd\ - 2r,|d|, sin0 + rjf , (B2) 

where dy, r y are given by 

rn = cos 0' cos + | sin 0' | y/l — cos 2 , (B3) 

d,l = rfcos0'. (B4) 

And finally, cos and sin 9 in terms of 9', 0' are given by 



cos (j) = d sin 2 0' + cos 0' y 1 — d? sin 2 0' , (B5) 

sin# = — (d n cos 2 9' + sin #'y/r 2 - d 2 cos 2 0> ) . (B6) 

The inverse transformation is given by interchanging 9 <-> 6 1 ', <-> 0' d <-> — d in the above 
expressions. 
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